/* Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010 Mario Orsi
   This file is part of Brahms.
   Brahms is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
   Brahms is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
   You should have received a copy of the GNU General Public License along with Brahms.  If not, see <http://www.gnu.org/licenses/>. */

/* Direct evaluation of diffusion [Rapaport, p. 124-126, 226].
   Sets of diffusion measurements are begun at fixed time intervals, so that several 
   sets of measurements based on different time origins will be in progress simultaneously
   because of the overlapped measurements.
   Memory is allocated in module startUp.c, function AllocArrays(). */

#include "dataStructs.h"

/* External data references: */
extern const VecR *siteCoords, *siteOrientVec; 
extern const Site *site; 
extern const VecR region; 
extern const int nSites, stepCount;
extern const real deltaT, timeNow;

// Global variables:
TBuf *tBuf; // buffer for time-dependent measurement [Rapaport, p. 510]
real *rrDiffuseAv; // translational diffusion mean [Rapaport, p. 128]
real *ddDiffuseAv; // rotational diffusion mean [Rapaport, p. 227]
int limitDiffuseAv; // total number of individual estimates used to produce the average value of the diffusion coefficient  
int nBuffDiffuse; // number of sets of data being collected at any time, each occupying a separate storage buffer 
int nValDiffuse; // number of measurements contributing to the set used to produce a single unaveraged estimate of the diffusion coefficient
int stepDiffuse; // measurement frequency (number of steps between measurements)

// Local (to this module) variables:
static int countDiffuseAv; // number of measurements in set [Rapaport, p. 510]

static void PrintDiffusion() // Adapted from Rapaport [p. 127]
{
  real measurementTime;
  int j;

  FILE *fPtr;
  char diffSnapFile[ 35 ];
  char time_ps[ 9 ]; // string containing timeNow, e.g. "100", in ps

  sprintf( time_ps, "%d", ( int ) ( 0.5 + timeNow ) );
  diffSnapFile[ 0 ] = '\0';
  strcat( diffSnapFile, "waterDiffusion_" );
  strcat( diffSnapFile, time_ps );
  strcat( diffSnapFile, "ps" );
  strcat( diffSnapFile, ".dat" );
  fPtr = fopen( diffSnapFile, "w" );
  for ( j = 0; j < nValDiffuse; j++ ) {
    measurementTime = j * stepDiffuse * deltaT;
    fprintf( fPtr, "%17.4f", measurementTime ); /* [ps], Brahms's fundamental unit of time (no need to convert) */
    fprintf( fPtr, "%12.5f", rrDiffuseAv[ j ] * Sqr( LENGTH_cm ) / TIME_s * 1e5 ); // [10^(-5)cm^2/s] 
    fprintf( fPtr, "%12.5f", ddDiffuseAv[ j ] ); 
    fprintf( fPtr, "\n" );
  }
  fclose( fPtr );
}

static void ZeroDiffusion() // [Rapaport, pp. 127, 227]
{
  int j;

  countDiffuseAv = 0;
  for ( j = 0; j < nValDiffuse; j++ ) {
    rrDiffuseAv[ j ] = 0.;  
    ddDiffuseAv[ j ] = 0.;
  }
}

/* InitDiffusion -- Initialises the calculation. Adapted from Rapaport [p. 127] 
   The values initially assigned to "count" determine the spacing between measurements: negative initial values delay
   the start of data collection for each set of measurements until the appropriate moment. */
void InitDiffusion()
{
  int nb;
  FILE *fPtr;

  fPtr = fopen( "brahms.log", "a" );
  fprintf( fPtr, "step %d -- initializing diffusion measurement\n", stepCount );
  fclose( fPtr );

  for ( nb = 0; nb < nBuffDiffuse; nb++ ) 
    tBuf[ nb ].count = - nb * nValDiffuse / nBuffDiffuse; 
  ZeroDiffusion();
}

static void AccumDiffusion() // [Rapaport, pp. 126-127, 226-227]
{
  real fac;
  int j, nb;

  for ( nb = 0; nb < nBuffDiffuse; nb++ ) {
    if ( tBuf[ nb ].count == nValDiffuse ) { /* if buffer "nb" is full */
      for ( j = 0; j < nValDiffuse; j++ ) { /* scan through the buffer elements */
	rrDiffuseAv[ j ] += tBuf[ nb ].rrDiffuse[ j ]; // accumulate translational diffusion
	ddDiffuseAv[ j ] += tBuf[ nb ].ddDiffuse[ j ]; // accumulate rotational diffusion
      }
      tBuf[ nb ].count = 0;  /* reset buffer counter */
      ++countDiffuseAv;
      if ( countDiffuseAv == limitDiffuseAv ) { /* if the #individual estimates required has been reached */

        fac = 1. / ( 6 * nSites * stepDiffuse * deltaT * limitDiffuseAv );
        for ( j = 1; j < nValDiffuse; j++ )
	  rrDiffuseAv[ j ] *= fac / j;

        fac = 1. / ( nSites * limitDiffuseAv );
        for ( j = 0; j < nValDiffuse; j++ )
           ddDiffuseAv[ j ] *= fac;

	PrintDiffusion();
	ZeroDiffusion();
      }
    }
  }
}

/* EvalDiffusion -- evaluate translational and rotational diffusion. Adapted from [Rapaport, p. 126, 226] */
void EvalDiffusion()
{
  VecR dr;
  int n, nb, ni;

  for ( nb = 0; nb < nBuffDiffuse; nb++ ) { /* scan through the individual buffers */
    if ( tBuf[ nb ].count == 0 ) { /* if empty buffer */
      DO_SITE {
        tBuf[ nb ].orgR[ n ] = siteCoords[ n ];
        tBuf[ nb ].rTrue[ n ] = siteCoords[ n ];
        tBuf[ nb ].orgD[ n ] = siteOrientVec[ n ];
      }
    }
    if ( tBuf[ nb ].count >= 0 ) {

      ni = tBuf[ nb ].count;

      tBuf[ nb ].rrDiffuse[ ni ] = 0.;
      DO_SITE {
        VSub( dr, tBuf[ nb ].rTrue[ n ], siteCoords[ n ] );
        VDiv( dr, dr, region );
        dr.x = Nint( dr.x );
        dr.y = Nint( dr.y );
        dr.z = Nint( dr.z );
        VMul( dr, dr, region );
        VAdd( tBuf[ nb ].rTrue[ n ], siteCoords[ n ], dr );
        VSub( dr, tBuf[ nb ].rTrue[ n ], tBuf[ nb ].orgR[ n ] );
        tBuf[ nb ].rrDiffuse[ ni ] += VLenSq( dr );
      }

      tBuf[ nb ].ddDiffuse[ ni ] = 0.;
      DO_SITE {
        tBuf[ nb ].ddDiffuse[ ni ] += VDot( tBuf[ nb ].orgD[ n ], siteOrientVec[ n ] );
      }
    }
    ++tBuf[ nb ].count;  /* increase buffer counter */
  }
  AccumDiffusion();
}

